% \input{gbmacros.tex}

\chapter{Inverse Problem}
\label{chap:inverse}

Cardiac motion estimation is the problem of determining a transformation that captures the motion of every point in the myocardium over the cardiac cycle. If we can register 3D images acquired at different phases of the cardiac cycle to each other, then we have an estimate of cardiac motion. The problem of motion estimation from cardiac cine images is ill posed, and relying solely on image similarity, even with very accurate similarity measures is not sufficient to capture the true motion of the heart. Current cardiac motion estimation methods rely on image similarity measure to drive the motion estimation, and typically incorporate a regularizer to smooth the deformation field. We propose to instead maximize the image similarity, subject to the motion estimate constrained by a mechanical model of the heart, which has been discussed in Chapter \ref{chap:model}.

The result of a 4D scan of a beating heart is a periodic sequence of $N$ 3-$D$ images, 

\[I(x,t)=\{ I_{t} (x),0\le t<N\} ,\] 

where $I_{0} $ is the end-diastolic image. We define the motion field as the transformation $\chi (x,t)$ defined over the image space. The transformation maps a point $x$ in the end-diastole frame $I_{0} $ of the image to its corresponding point in frame $I_{t} $ at time $t$. This is illustrated in Figure \ref{fig:motionForm}. The displacement field $U(x,t)$ defines the mapping from the coordinate system of the end-diastole image $I_{0} $ to the image at time $t$, $I_{t} $. The transformation and the displacement are related as, $\chi (x,t)=x+U(x,t)$. 

\begin{figure}
\begin{center}
	\includegraphics[width=\textwidth]{images/motion-est.pdf}
 \end{center}
 \caption{The formulation of the motion estimation problem.}
 \label{fig:motionForm}
\end{figure}

\section {Formulation of the inverse problem}
\label{sec:model}

The basic premise of our formulation is the following: The heart
motion is induced by the active forces in the myocardium. If we knew
the exact biomechanical model for the myocardial tissue (constitutive
law, geometry, fiber orientations, material properties for the heart
and surrounding tissues, endocardial tractions due to blood flow) and
the active stretching time-space profile, then we could solve the
so-called ``{\em forward problem}'' for the displacements of the
myocardial tissue. Similarly, if we knew the displacements at certain
locations in the myocardium, we could solve the so-called ``{\em
inverse problem}'' to reconstruct active forces so that the motion due
to the reconstructed forces matches the observed one. More generally,
we have cine-MRI data but not the displacements. We can still invert
for the displacements---by solving a biomechanically-constrained image
registration problem.

In this context, an abstract formulation of the myocardium motion
estimation problem is given by
%%
\begin{equation}
\label{e:a} % e:abstract
  \min_{u,s}  \MA{J}(I_t,I_0,u)\quad   \mbox{subject to}\quad \MA{C}(u,s) = 0.  
\end{equation}
%%
Here, $I_t:=I_t(x,t)$ is the cine-MR image sequence with $x,t$
denoting the space-time coordinates, $I_0:=I(x,0)$ is the initial
frame (typically end-diastole), $u:=u(x,t)$ is the {\em displacement}
(motion), $s=s(x,t)$ is the active fiber contraction, and $\MA{C}$ is
the forward problem operator. Also, $\MA{J}$ is an image similarity
measure functional.  This is a classical PDE-constrained inverse
problem \cite{gordon03}. Notice that there is no need
for elastic, fluid, or any kind of regularization for $u$. It is
constrained through the biomechanical model
$\MA{C}$.\footnote{However, one can show that the problem is ill-posed
on $s$. Here we regularize by discretization of $s$.}


\subsection{Objective function ($\MA{J}$)} 

For the purposes of this thesis, we assume a preprocessing step, in which the point correspondences are computed manually for all frames, i.e.,
${d_j(t):=u(x_j,t)}_{i=1}^M$ at $M$ points. Then, the objective function is given by
%%
% \int_{#2}^{#3} #4 \, d\: #1}
\begin{equation} 
\label{e:o}
 \MA{J}:=\ivar{t}{0}{1}{(Qu-d)^2}:=\ivar{t}{0}{1}{\sum_{i=1}^M (u(x_j,t)-d_j(t))^2},
\end{equation}

where $Q$ is the so called spatial observation operator. 


\subsection{Inverse problem} 
The inverse problem is stated by \eqref{e:a}
where $\MA{J}$ is given by \eqref{e:o} and $\MA{C}$ is given by
\eqref{e:fiberForce}. By introducing Lagrange multipliers $p$, the first-order
optimality conditions for $\eqref{e:a}$ can be written as:
%%
\begin{equation} 
\label{e:kkt}
\begin{split}
M\ddot{u}(t) + C\dot{u}(t) + K u(t) + A s(t) = 0, \quad \dot{u}(0)
=u(0) = 0,\\ 
M\ddot{p}(t) - C\dot{p}(t) + K p(t) + Q^T(Qu-d) = 0, \quad \dot{p}(1) =
p(1) = 0,\\ 
A^T p(t) = 0.  
\end{split}
\end{equation}
 
%%
The second equation is the so-called ``{\em adjoint
problem}''. Equation \eqref{e:kkt} consists of a system of
partial-differential equations for $u$ (cardiac motion), $p$
(adjoints), and $s$ (active fiber contraction). It is a 4D boundary
value problem since we have conditions prescribed at both $t=0$ and
$t=1$. 


%\subsection{Discretization and solution algorithms}
\subsection{Discretization and solution algorithms}  
We discretize the
forward and adjoint problems in space using a Ritz-Galerkin
formulation. We have developed a parallel data-structure and meshing
scheme, discussed in \cite{sundar-sampath-biros-e07}. The basis
functions are trilinear, piecewise continuous polynomials. In time, we
discretize using a Newmark scheme. The overall method is second-order
accurate in space and time.  The implicit steps in the Newmark scheme
are performed using Conjugate Gradients combined with a
domain-decomposition preconditioner in which the local preconditioners
are incomplete factorizations. The solver and the preconditioner are 
part of the PETSc package \cite{petsc-web-page}.

For these particular choices of objective function and forward problem
the inverse problem \eqref{e:kkt} is linear in $p$, $u$, and $s$. We
use a reduced space approach in which we employ a matrix-free
Conjugate-Gradients algorithm for the Schur-complement of $s$---also
called the (reduced) Hessian operator. Each matrix-vector
multiplication with the Hessian requires one forward and one adjoint
cardiac cycle simulation. Furthermore, one can show that the Hessian
is ill-conditioned. Thus, the overall computational cost of the
inversion is high.  We are developing efficient multigrid schemes for
the inverse problem. Details of this approach can be found in
\cite{gordon03}. To reduce the computational cost for
the calculations, we used a reduced-order model
for $s$ in which $\psi$ is a product of B-splines in time and radial
functions in space (Gaussians).  This discretization not only does it
allow acceleration of the inverse problem but it introduces a model
error since the synthetic ``ground truth'' is generated using a full
resolution fiber model and the inversion is done using the reduced
resolution fiber model. This allows to perform preliminary tests on
the sensitivity of our method to model errors.
%}}}


\subsection{Active Force Models}

We consider two models for the active forces: the first being the generic case where arbitrary forces can be specified and the second model which accounts for the contractility of the myocytes. For the generic case the force term is given by,

\begin{equation}
F(t) = M\grbf{f}(t),
\end{equation}

where, $\grbf{f}(t)$ are the nodal forces at time $t$.

To model the contractility of the myocytes given the fiber contractility $\grbf{s}$ as a function of space and time and the myocyte orientation $n$, we define the active stretch tensor $U=1+ \grbf{s}\, n\otimes n$, whose divergence results in a distributed active force of the form $\Div(\grbf{s}\, n\otimes n)$. The force term in this case is given by, 

\begin{equation}
 F(t) = A\grbf{s}(t), \quad A_{ij}= \int (n \otimes n) \Grad \phi_i \phi_j.
\end{equation}
 
To reduce the computational cost for the calculations, we used a reduced-order model for $\grbf{s}$ and $\grbf{f}$ as a combination of B-splines bases in time and radial basis functions in space (Gaussians). The forces can then be written in terms of the B-spline basis, $B$ and the radial basis, $G$, and the control parameters $\grbf{\mu}$,

\begin{eqnarray}
	\grbf{f}(\grbf{x}, t) &=& \sum_{k=1}^3 \grbf{e}_k \sum_i G_i^k(\grbf{x}) \sum_j \grbf{\mu}_{ij} B_{j}(t), \\
	\grbf{s}(\grbf{x}, t) &=& \sum_i G_i(\grbf{x}) \sum_j \grbf{\mu}_{ij} B_{j}(t).
\end{eqnarray}

In the matrix form this defines the parametrization matrix $C$ is given by,

\begin{equation}
\label{e:C1}
C_{xt, ij} = G_i(x) B_{j}(t).
\end{equation}

We can write the active forces in terms of the parametrization matrix $C$, as

\begin{eqnarray*} 
F &=& M \grbf{f} = M C \grbf{\mu}, \\
F &=& A \grbf{s} = A C \grbf{\mu}.
\end{eqnarray*}


Additionally we need to define the transpose of $A$ and $C$ since they appear in the reduced gradient and reduced Hessian operators. The transpose of the myocyte parametrization matrix $A$ is given by,

\begin{equation}
\label{e:Atrans}
A^T_{ij}= \int \Grad \phi_j \phi_i (n \otimes n).
\end{equation}
 
The transpose of the parametrization matrix $C$ is given by,

\begin{equation}
C^T_{ij, xt} = G_i B_{j}
\end{equation}

%}}}


%{{{ Results 

\section{Results}

In this section we describe experiments conducted to validate the motion estimation algorithm. We conducted experiments using synthetic models to determine the numerical correctness and robustness of the estimation with respect to noise and the estimation parameters. These are described in Section \ref{sec:synResults}. We also tested the estimation using tagged MR datasets to validate the effectiveness of the method on real data. These experiments are described in Section \ref{sec:taggedResults}. 

\subsection{Synthetic Datasets}
\label{sec:synResults}

In order to assess the fiber-orientation parametrized model of the forces (\ref{e:fiberForce}), we use an ellipsoidal model of the left ventricle. The fiber orientations are generated by varying the elevation angle between the fiber and the short axis plane between $+60^{\circ}$ and $-60^{\circ}$ from the endocardium to the epicardium \cite{sachse04,sermesant04}. The model along with the fiber orientations is shown in Figure \ref{f:fiberModel}. For this model we selected a Poisson's ratio $\nu=0.45$ and a Young's modulus of 10 kPa for the myocardial tissue and 1 kPa for the surrounding tissue and ventricular cavity. Raleigh damping ($C =\alpha M + \beta K$) was used with parameters $\alpha=0$ and $beta=7.5\times10^{-4}$. In order to drive the forward model, we generated
forces by propagating a synthetic activation wave from the apex to the
base of the ventricles. Snapshots of this activation wave at different
phases of the cardiac cycle are shown in Figure \ref{fig:activation}.


\begin{table}
  \centering
  \begin{tabular}{|c|c|c|}
		\hline
		Basis Size & Relative Error $\|\cdot\|_2$ & Time \\
		\hline
		$2^3$ &  $1.31\times10^{-1}$ & 36 mins \\
		$4^3$ &  $5.67\times10^{-2}$ & $\approx$ 5 hrs  \\
		$4^3$ &  $1.12\times10^{-1}$ & 108 mins  \\
		$8^3$ &  $9.66\times10^{-2}$ & 141 mins \\ 
		\hline	
		\end{tabular}
	\caption{Error in recovery of activation for
  increasing number of radial basis functions. By changing the
  inversion solver accuracy, we can accelerate the calculation without
  compromising accuracy (e.g., the $4^3$ calculation). }
	\label{tab:invRecovery}
 \end{table}

\begin{table}
  \centering
  \begin{tabular}{|c|c|}
		\hline
		Observations & Relative Error $\|\cdot\|_2$ \\
		\hline
		Full &  $5.36\times10^{-2}$  \\
		$1/8$ &  $6.21\times10^{-2}$  \\
		$1/64$ &  $8.51\times 10^{-2}$ \\ 
		\hline	
		\end{tabular}
	\caption{Error in the recovery of activation with partial
	observations of the displacements. Errors are reported on the
	cylinder model for a grid size of $32$ with $4^3$ basis functions.}
	\label{tab:invPartial}
\end{table}


For the inverse problem, we validate the error in the estimation of
the activations for different degrees of parametrization using the
radial basis. In all cases, the number of b-Spline basis per spatial
location were fixed to 5 degrees of freedom. The relative error in the
estimation of the activation for a $64^3$ grid is for spatial
parametrizations of $2^3$, $4^3$ and $8^3$ is tabulated in Table
\ref{tab:invRecovery}.  Ground truth activations are generated using a wave propagating from the apex to the base of the ventricle \ref{fig:activation}. These runs were done on 64 processors.  In
addition, we investigated the error in the estimation when only
partial observations are available. We compared estimations based on
full and sparse observations with $12\%$ and $6\%$ samples against the
analytical solutions. These results are tabulated in Table
\ref{tab:invPartial}.
%% if time left...
In order to assess the sensitivity of the motion estimation framework,
we estimated the motion for the synthetic model of the heart at a grid
size of $64$ with a radial basis parametrization of $4^3$ by adding
noise to the system. We added a $5\%$ random error on the estimates of
the fiber orientation and to the material properties of the
myocardium. In addition, we added a $1\%$ noise to the true
displacements. The system converged and the relative error ($L_2$)
increased from $5.67\times10^{-2}$ to $9.43\times10^{-2}$.

\subsection{Validation using tagged MR images}
\label{sec:taggedResults}

\subsubsection{Data}

We acquired cine MR sequences with and without myocardial tagging for 5 healthy volunteers in order to validate our motion estimation algorithm. The cine MR and tagged images were acquired on a Siemens Sonata 1.5T\texttrademark~  scanner during the same scanning session at end-expiration, thus the datasets are assumed self-registered. Short and long axis segmented k-space breath-hold cine TrueFISP (SSFP) images were acquired for three healthy volunteers. The image dimensions were 156x192 pixels with a pixel size of 1.14x1.14 mm2. A total of 11 slices were acquired with a slice thickness of 6mm and a slice gap of 2mm. For validation purposes, we also acquired grid tagged, segmented k-space breath-hold cine TurboFLASH images. The image dimensions were 156x192 pixels with a pixel size of 1.14x1.14 mm2. The slice thickness was 8mm. Three short axis (apical, mid-ventricular and basal) and one long-axis images were acquired. The displacements (observations) at the tag intersection points within the myocardium were computed manually using the markTags GUI which allows users to manually place points at the tag intersections. The program also allows the users to define the correspondence over time and is capable of generating temporally interpolated displacements fields needed for the motion estimations. The displacements were interpolated in time using linear interpolation. A screenshot of the markTags program is shown in Figure \ref{fig:markTags}. An average of 70 tag intersection points were selected over the left and right ventricles on each plane resulting in around 300 observations in space. Three independent observers processed all 5 datasets to get three sets of observations. The average displacements were used as the true observations for the motion estimation algorithm.

\begin{figure} 
\centering
\includegraphics[width=0.6\textwidth]{images/markTags}
\caption{The {\tt markTags} program that is used to manually track the tag intersection points.}
\label{fig:markTags}
\end{figure}

%}}}

\subsubsection{Pre-processing}

The patient specific fiber orientations are required for the motion estimation algorithm. The end-diastolic MR image was segmented manually using the Insight SNAP application. The segmented end-diastolic image was registered to the MR image of an ex-vivo human heart (atlas) for which diffusion tensor imaging was performed. We used the in-house HAMMER algorithm \cite{hammer} for performing the non-rigid registration with the atlas and template warping was used to estimate the patient specific myocardial fiber orientation as explained in Section \ref{sec:warpDTI}. The segmentations were also used to assign material properties to the myocardial tissue, the blood and the surrounding tissue.

\subsubsection{Experiments}

For each of the 5 subjects, we had observations from 3 independent observers in the form of manually placed landmarks over the entire image sequences. To account for observer variability, we treated the mean location of the landmark as the ground truth observations. We used a grid size of $64^3$ for all inversions. In three seperate experiments we used $70\%$, $50\%$, and $20\%$ of the ground truth observations as the data for the inversion. The observations that are not used during the inversion are the control observations and are used for assessing the goodness of the inversion.

We used the fiber orientation parametrized force model \ref{e:fiberForce} along with the B-spline and the radial bases to reduce the paramter space. We used a total of $4^3$ spatial parameters and for each 5 B-Spline temporal parameters, giving us a total of 320 parameters.     

After inversion, one additional forward solve was performed using the dense estimates of the fiber contractions to obtain dense estimates of myocardial displacements. These displacements were compared with the ground truth displacements at the tag-intersection points. The relative error (as a percentage) and the absolute error in millimetres for all observations and restricted to only the control observations are shown in Table \ref{tab:tagged}. A visual comparison of the estimation with different levels of partial observations overlaid on the original tagged images are shown in Figure \ref{f:tagged}.


\begin{table}
  \centering
  \begin{tabular}{|c|c|c|c|c|}
		\hline
		 & Error ($\%$) & Error ($\%$) &  Error ($mm$) &  Error ($mm$) \\
		 & All obs. & Control obs. & All obs. & Control Obs. \\   	
		\hline
		 Observers & 9.21 & 8.41 &	1.39 & 1.26\\
		 Algorithm $70\%$ & 12.46 &14.22 & 1.88 & 2.13 \\
		 Algorithm $50\%$ & 16.37 & 21.02 & 2.47 & 3.15 \\ 
		 Algorithm $20\%$ & 41.61 & 51.19 & 6.28 & 7.67 \\
		\hline	
		\end{tabular}
	\caption{Results from tagged data }
	\label{tab:tagged}
 \end{table}

\begin {figure}[tbp]
\centering
	\subfloat[Manual]{
		\includegraphics[width=0.45\textwidth]{images/estimate_manual}
	}
	\subfloat[Estimate with $70\%$]{
		\includegraphics[width=0.45\textwidth]{images/estimate_70}
	} \\
	\subfloat[Estimate with $50\%$]{
		\includegraphics[width=0.45\textwidth]{images/estimate_50}
	} 
	\subfloat[Estimate with $20\%$]{
		\includegraphics[width=0.45\textwidth]{images/estimate_20}
	}

\caption{Comparison of the motion estimation algorithm using partial observations.}
\label{f:tagged}
\end{figure}

\begin {figure}[tbp]
\centering
	\subfloat[end-diastole] {
		\includegraphics[width=0.3\textwidth]{images/def01} }
		\subfloat[mid-systole] {
		\includegraphics[width=0.3\textwidth]{images/def02} }
		\subfloat[end-systole] {
		\includegraphics[width=0.3\textwidth]{images/def03} }
\caption{The deformation field shown for a mid ventricular slice for (a) end-diastole, (b) mid-systole, and (c) end-systole.}
\label{f:heartDefs}
\end{figure}

The deformation fields at end-diastole, mid-systole and end-systole for the mid ventricular slice from one subject is shown in Figure \ref{f:heartDefs}.  

An alternate way of interpreting the results is to analyze the left-ventricular volume, obtained using manual segmentation and by the using the motion estimates\footnote{using $70\%$ observations}, as seen for one particular subject in Figure \ref{f:heartVols}. The motion estimation based volume curve is smoother than the one obtained by manual segmentation, but is in general agreement with the latter. We manually segmented the left ventricle for the 5 subjects and evaluated the difference in ventricular volume over the entire cardiac cycle. The average volume error was $3.37\%$, with standard deviation $2.56\%$; and average volume overlap error was $7.04\%$, with standard deviation $3.28\%$.

\begin {figure}[tbp]
\centering
	\includegraphics[width=0.8\textwidth]{images/heartVols}
\caption{Left ventricular volume of a selected subject, segmented by our algorithm (solid curve) and by hand (dotted curve) over all frames in a cardiac cycle.}
\label{f:heartVols}
\end{figure}


 \section{Conclusions}
We presented a method for cardiac motion reconstruction. We integrate
cine-MR images and a biomechanical model that accounts for
inhomogeneous tissue properties, fiber information, and active
forces. We presented an inversion algorithm. We will able to solve
problems that involve 300 million unknowns (forward and adjoint in
space-time) in a couple of hours on 64 processors---a modest computing
resource. The potential of the method as multicore platforms become
mainstream is significant.

The limitations of our current implementation (but not the method) is
the assumptions of linear geometric and material response and the
potential bias due to template-based fibers that does not account for
anatomical variability, that is still requires some preprocessing of
the initial frame to assign material properties and fiber orientation,
that assumes zero residual stresses and initial conditions, and that
it does not include an electrophysiology model.

Our on-going work includes transition to an intensity-based
image-registration inversion (in which case we need to solve a
nonlinear inversion), further verification of the method by comparing
to manually-processed real data, and its clinical validation by
reconstructing motions of normal and abnormal populations and
conducting statistical analysis.  Among the many open problems are the
level of required model complexity for clinically relevant motion
reconstructions, the bias of the fibers, the sensitivity to the values
of the material properties, and the sensitivity to the image
similarity functional.





